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HIGHLIGHTS 


•  A  radial  mass-transport  coefficient  is  used  to  modify  the  thermal-electrochemical  Dualfoil  model. 

•  Simulation  of  current  and  temperature  under  very-high-current  discharge  is  achieved. 

•  Mass-transport  limitation  decreased  by  high  cell  temperature  causes  a  rapid  temperature  rise. 

•  Effective  heat  transfer  outside  the  cell  center  is  critical. 
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The  Dualfoil  model  is  used  to  simulate  the  electrochemical  behavior  and  temperature  rise  for  MCMB/ 
LiCo02  Li-ion  cells  under  a  small  constant-resistance  load,  approaching  a  short-circuit  condition.  Radial 
mass  transport  of  lithium  from  the  center  of  the  pore  to  the  pore  wall  has  been  added  to  the  model  to 
describe  better  current  limitations  at  very  high  discharge  currents.  Electrolyte  and  solid-surface- 
concentration  profiles  of  lithium  ions  across  the  cell  at  various  times  are  developed  and  analyzed  to 
explain  the  lithium-ion  transport  limitations.  Sensitivity  tests  are  conducted  by  changing  solution  and 
solid-state  diffusion  coefficients,  and  the  heat-transfer  coefficient.  Because  diffusion  coefficients  increase 
at  high  temperature,  calculated  discharge  curves  can  show  currents  dropping  initially  but  then  rising  to  a 
second  peak,  with  most  of  the  available  capacity  being  consumed  in  the  second  peak.  Conditions  which 
lead  to  such  a  second  peak  are  explored. 

Published  by  Elsevier  B.V. 


1.  Introduction 

Li-ion  batteries  are  widely  used  in  mobile  technologies  because 
of  their  high  energy  and  power  density.  The  usage  extends  from 
small  cell  phones  to  electric  vehicles,  and  even  to  airplanes.  With 
bigger  batteries,  more  concern  on  battery  safety  is  aroused.  Re¬ 
views  of  thermal  issues  in  Li-ion  batteries  point  out  that  the  safety 
problem  results  from  thermal  runaway  causing  fire  and  explosion, 
and  it  is  clear  that  no  existing  thermal-management  strategy  or 
technology  meets  all  these  requirements  [1,2], 

The  key  to  battery  safety  is  to  manage  the  temperature  properly. 
A  battery  can  be  safe  under  moderate  operating  condition  (e.g.,  low 
current),  but  is  compromised  for  larger  cells  and  higher  currents 
(low-resistance  loads).  The  temperature  rise  of  a  Li-ion  cell  is 
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determined  by  heat  produced  internally  and  heat  exchanging  with 
the  environment  during  operation,  and  can  be  predicted  using  a 
thermal-electrochemical-coupled  model.  The  coupling  of  thermal 
and  electrochemical  models  arises  from  temperature-dependent 
physicochemical  properties  and  the  heat-generation  equation. 
The  history  of  thermal-electrochemical-coupled  models  has  been 
given  in  Refs.  [3,4],  Internal  and  external  short-circuit  tests  are 
routine  for  Li-ion-battery  safety  tests  in  industry,  and  high  tem¬ 
peratures  that  result  from  extremely  high  discharge  current  are 
sometimes  observed  during  short-circuit  tests  [5—8],  It  is  urgent  to 
interpret  the  various  experimental  phenomena  fundamentally  by 
simulation.  The  goal  of  this  work  is  to  use  the  thermal- 
electrochemical-coupled  Dualfoil  model  to  simulate  the  tempera¬ 
ture  rise  in  Li-ion  cells  at  very  high  currents  and  understand 
fundamentally  the  quick  temperature  rise  to  above  120  °C  during 
short-circuit  tests. 

The  Dualfoil  model  is  a  powerful  and  flexible  mathematical 
model  that  can  be  used  to  treat  the  coupled  phenomena  in  a 
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porous-electrode  battery  system  (including  Ni-MH  and  Li-ion 
batteries).  It  has  six  coupled  differential  equations  (a  material 
balance  on  the  electrolyte,  a  material  balance  on  the  solid  inter- 
calant,  Ohm's  law  in  the  liquid  phase,  Ohm's  law  in  the  solid  phase, 
a  current  balance  that  relates  the  flow  of  current  between  the  solid 
and  liquid  phases,  and  a  Butler— Volmer  kinetics  expression)  [9], 
Review  of  the  Dualfoil  model  and  more  detailed  descriptions  are 
given  in  Refs.  [10—12],  A  clear  schematic  of  the  Li-ion  cell  in 
Dualfoil  modeling  has  been  given  in  a  recent  reference  [13].  Unlike 
other  commercial  simulation  programs,  such  as  COMSOL,  MATLAB, 
the  Dualfoil  program  based  on  Dualfoil  model  is  open  access,  and 
can  be  downloaded  from  our  website  and  has  low  computational 
demands  [14].  The  Dualfoil  software  employs  Newman’s  BAND  (j) 
subroutine,  which  is  a  finite-difference  method  (FDM  )  used  to 
simulate  electrochemical  systems  for  more  than  four  decades.  The 
model  has  been  improved  over  the  years  to  include  impedance 
[15,16],  a  better  thermal  balance  [17—21],  additional  battery 
chemistries  [22—25],  etc.,  and  has  more  and  more  applications 
[13,26—29],  When  discharging  under  nonequilibrium  conditions 
such  as  very  high  currents,  the  distribution  of  current,  potential, 
and  lithium  concentration  inside  the  porous  electrode  is  nonuni¬ 
form,  leading  to  increased  internal  resistance,  producing  large 
heat,  and  thus  increasing  the  cell  temperature.  The  transport 
properties  in  the  electrolyte  and  solid  phases  and  the  electro¬ 
chemical  reaction  rate  can  be  increased  by  temperature  rise 
following  some  empirical  relationship,  and  an  Arrhenius  expo¬ 
nential  relationship  is  usually  adopted.  This  temperature  rise  in 
turn  will  cause  more  current  to  pass  the  cell.  This  kind  of  positive 
feedback  between  temperature  and  those  physicochemical  prop¬ 
erties  can  lead  to  overheating  issues  and  must  be  dealt  with  30], 
Simulation  can  give  deep  insight  into  the  internal  behavior  of 
batteries,  and  can  provide  information  that  is  difficult  or  impos¬ 
sible  to  obtain  by  experiment  [31],  A  good  review  of  model 
development  and  simulation  of  Li-ion  cells  at  different  length  and 
time  scales  is  given  by  reference  [29], 

When  simulating  the  thermal  and  electrochemical  behavior  of 
a  Li-ion  cell  at  abnormally  large  discharge  currents,  there  are  two 
difficult  problems  for  the  Dualfoil  model.  The  first  problem  is  that 
the  salt  concentration  can  be  reduced  substantially  and  even  try 
to  go  negative  locally  when  the  discharge  current  is  very  high, 
and  the  Dualfoil  program  has  shown  a  tendency  to  stop  short  of 
the  full  simulation  due  to  lack  of  convergence.  For  a  real  elec¬ 
trolyte  system,  the  local  salt  concentration  cannot  be  zero  or 
negative.  This  issue  has  always  been  present  in  the  model,  and 
various  approximations  have  been  developed  to  bypass  the  issue 
for  a  given  set  of  conditions.  In  this  work,  we  develop  a  radial 
mass-transport  coefficient  of  solution  diffusion  to  solve  the 
convergence  problem.  Most  previous  porous-electrode  models 
consider  only  mass  transport  along  the  pores  and  ignore  radial 
transport  in  the  pores.  This  means  that  the  absence  of  a  radial 
transport  limitation  allows  even  the  most  dilute  solutions  to 
support  large  currents.  To  address  this  issue,  a  radial  mass- 
transport  condition  (i.e„  limiting  current)  has  been  introduced 
into  the  model.  With  discharge  at  extremely  high  currents,  the 
local  salt  concentration  can  only  get  very  close  to  zero,  and  the 
transfer  current  will  approach  its  limit.  This  modification  to  the 
Dualfoil  model  results  in  smooth  operation  and  allows  the  pro¬ 
gram  to  run  at  virtually  all  currents,  and  dramatically  reduces  the 
running  time.  The  model  modification  and  derivation  are  pre¬ 
sented  in  the  next  section. 

The  second  problem  is  how  to  compute  the  variable  solid 
diffusion  for  lithium  ions  with  more  accuracy  and  in  short 
computational  time.  In  this  work,  lithium  ions  are  still  assumed  to 
diffuse  radially  in  the  spherical  solid  particle  following  Fields  laws; 
and  with  a  solid  diffusivity  which  does  not  vary,  the  concentration 


profile  can  be  computed  efficiently  by  the  Duhamel  superposition 
integral  (which  is  an  advanced  mathematical  method  [32])  as  done 
in  the  original  Dualfoil  model.  There  are  many  mathematical 
methods  for  treating  solid  diffusion  in  simulation  research,  and  the 
usage  of  some  approximate  approaches  for  simulations  of  Li-ion 
solid  transport  in  batteries  is  reviewed  133],  More  detailed  dis¬ 
cussion  of  the  advantages  and  disadvantages  of  using  Duhamel 
superposition  integral  to  compute  with  a  variable  solid  diffusion  is 
shown  in  a  later  section. 

In  this  paper,  the  thermal-electrochemical  model  (Dualfoil)  is 
used  to  simulate  the  electrochemical  behavior  and  temperature  rise 
for  an  MCMB/LiCo02  cell  discharging  under  a  constant-resistance 
load.  A  set  of  parameters  is  used  in  the  base-case  simulation  for  a 
small  resistance  load  (Rioad)  of  0.2  x  1CT3  0  m2.  Then  sensitivity 
tests  on  some  parameters  are  conducted  and  investigated.  By 
looking  at  very  high  currents  (approaching  a  short-circuit  condi¬ 
tion),  we  hope  to  learn  more  about  parameters  which  can  impact 
overheating  in  batteries. 


2.  Model 


In  the  Dualfoil  model,  the  six  coupled  differential  equations  are 
solved  numerically  and  simultaneously  at  each  time  step.  A  sum¬ 
mary  of  the  equations  is  presented  in  Table  5  in  the  appendix;  it  is 
noted  that  only  the  kinetic  equation  is  modified  by  a  denominator. 
A  function  has  been  developed  which  describes  the  radial  mass 
transport  of  lithium  ions  in  the  pore.  Starting  with  the  original  ki¬ 
netic  B— V  (Bulter— Volmer)  equation  in  the  model,  (i  is  transfer 
current  density,  j  is  transfer  flux) 


1  =  jF  =  i'o 


exp 


acF  \] 
RT  Vs) 


[1] 


insert  a  new  term  in  front  of  the  cathodic  term 
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Physically  this  is  intended  to  represent  a  term  like  co/cbuik  =  1-1/ 
him,  where  ium  is  proportional  to  the  bulk  or  average  concentration 
in  the  pore  and  can  be  calculated  roughly  by  equation  (3).  c  is  the 
average  concentration  already  calculated  by  Dualfoil  in  the  pores, 
and  r  is  the  radius  of  the  pore,  which  represents  a  diffusion-layer 
thickness  and  can  only  be  estimated.  At  the  cathode,  the  transfer 
current  is  negative,  and  hjm  should  be  negative. 
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I  him  I  —  r 


[3] 


The  new  term  denominator  is  added  in  the  new  B— V  equation  as 
following. 
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and  coef  represents  -io/Oiim/c)-  As  io,  i[im  are  all  concentration 
dependent,  and  iiim  can  only  be  estimated,  we  set  coef  as  constant  in 
the  program  to  simplify  this  problem. 

This  is  not  intended  to  be  a  rigorous  addition  to  the  program; 
it  is  designed  to  make  convergence  much  easier  while  still  being 
based  on  a  physical  concept.  The  exact  value  assigned  to  coef 
does  not  have  a  great  significance  in  terms  of  the  results.  This 
idea  provides  a  physically  realistic  modification  to  prevent 
driving  the  solution  concentration  locally  to  zero.  Then  the 
Dualfoil  program  is  able  to  run  smoothly  and  quickly  under  very- 
high-current  discharges  (in  excess  of  13000  A  rrT2).  As  the  so¬ 
lution  concentration  approaches  zero  in  the  pore,  the  mass 
transport  is  slowed  down  due  to  limited  rates  of  radial  diffusion. 
The  mass  transfer  is  not  infinite  but  approaches  a  limit,  and  the 
transfer  current  i  approaches  the  limiting  current  (fiim).  With 
r  =  5  x  1CT6  m,  c  =  0.01  mol  rrT3,  and  D  =  2.0  x  1CT10  m2  s_1,  iiim 
is  calculated  as  0.04  A  m~2.  By  changing  coef  in  the  model,  the 
transfer  current  in  the  back  of  positive  electrode  should  be 
approaching  iiim.  Intensive  simulation  results  indicate  that 
coef  <  1  mol  m  3  can  achieve  the  goal.  In  the  following  simula¬ 
tions,  coef  is  set  to  1  mol  m  3.  To  see  more  details  on  this 
modification  (denominator,  coef),  one  can  explore  the  latest 
version  of  Dualfoil  program  on  our  website  [14], 

3.  Model  parameters 

For  the  model  based  on  fundamental  laws  of  transport,  kinetics, 
and  thermodynamics,  many  physical  properties,  as  well  as  the  ge¬ 
ometry  of  the  cell,  are  required.  The  parameters  used  in  the  base 
case  are  listed  in  Tables  1—4.  The  limit  on  number  of  iterations  is  set 
to  a  relatively  large  number  so  that  convergence  can  be  achieved  on 
a  zero-time  step.  The  numbers  of  nodes  in  the  negative  electrode, 
separator,  positive  electrode,  and  solid  particle  are  set  to  be  40,  40, 
80,  and  100.  The  discharge  cutoff  potential  is  0.0001  V.  The  resis¬ 
tance  load  is  very  small,  that  is,  Rioad  =  0.0002  fi  m2. 

In  this  paper,  we  use  the  simple  one-dimensional  (ID)  thermal 
and  electrochemical  model  for  one  sandwich  cell,  the  detail  of  the 
heat  conduction  within  the  cell  are  not  discussed,  and  the  tem¬ 
perature  distribution  in  the  cell  is  assumed  to  be  uniform.  Only  the 
heat  dissipated  to  the  surroundings  from  the  cell  surface  is 
considered.  This  simplification  is  realistic  and  reasonable  for  rela¬ 
tively  small  cells  of  elongated  geometry  which  have  a  high  aspect 
ratio,  such  as  a  prismatic  cell  with  stacked  layers,  or  a  cylindrical 
cell  having  spirally  wound  layers. 


Table  1 

Design-adjustable  parameters. 


Parameters 

UC6  (MCMB) 

LiyCo02 

Thickness  of  electrode  (pm) 

96 

60 

Thickness  of  current  collector  (pm) 

10 

10 

Initial  stoichiometric  value  (xQ  or  yQ) 

0.8 

0.6 

Particle  radius  (pm) 

8 

5 

Volume  fraction  of  electrolyte,  e 

0.4 

0.36 

Volume  fraction  of  inert  filler,  ef 

0.064 

0.106 

Parameters 

Value 

Ambient  air  temperature  (K),  Ta 

298 

Volume  fraction  of  electrolyte  in  separator 

0.4 

Separator  thickness  (pm) 

25 

Grid  resistance  (Q  m2),  RG 

0.0008' 

Grid  resistance  outside  the  cell  sandwich  (Q  m2),  RGext 

0.00023 

Mass  of  the  whole  cell  (kg  m~2),  M 

0.4932 

Heat  capacity  of  system  (J  kg-1  K),  Cp 

1000b 

Heat-transfer  coefficient  (W  m-2  K),  h 

0.368L 

a  ref.  [36], 
b  ref.  [4,34], 
c  ref.  [35], 


Table  2 

Physical  property  values  at  298  K. 


Parameters 

LixCg  (MCMB) 

LiyCo02 

Matrix  conductivity  (S/m),  o 

1003 

0.53 

Film  resistance  (Q*m2),  Rf 

0.0035b 

0 

Coulombic  capacity  (mAh/g) 

372 

148 

Maximum  Li+  concentration  (mol  m-3),  ct 

2.86  x  10“ 

2.55  x  104 

Density  of  active  material  (kg  m~3) 

1 800c 

501 0C 

Diffusion  coefficient  in  solid  (m2  s-1),  Ds 

7  x  1(T14cI 

3.0  x  10-14d 

Transfer  coefficient  aa,  ac 

0.5 

0.5 

Reaction  rate  constant  (m2  5  mol-0  5  s-1),  k 

3.0  x  l(T9a 

3.0  x  l(T9a 

Density  of  current  collector  (kg  m-3) 

8954c 

2707c 

a  Assumed  to  be  large. 
b  Ref.  [24], 
c  Ref.  [37], 
d  Assumed. 


Table  3 

Parameters  for  the  electrolyte  (1  M  LiPF6  in  1:1  EC/DMC).a 


Parameters  Value 


Density  (kg  m-3) 

1324 

Salt  concentration 

1000 

(mol  m-3),  c 

Salt  diffusion  coefficient 

5.34  x  10  10E15exp( 

(m2  s-1),  D 

Lithium  ion  transference 

0.4 

number:  tj 

Ionic  conductivity  of 
liquid  solution  (S/m),  k 

^5[00911  +  1i9o1oT-1  052(iooo) 

+  o.l  554 

Viooo/  J 

Thermodynamic  factor: 

1 

(1  +  d  In  f±/d  In  c) 

3  use  the  property  put  in  the  Dualfoil  program  directly. 


The  heat  capacity,  Cp,  is  estimated  to  be  1000  J  kg-1  -K,  including 
all  the  cell  materials  [4,34].  The  capacity  ratio  of  negative/positive 
electrode  is  1.56.  The  cell  is  positive  limited,  and  its  capacity  is 
calculated  to  be  17.6  Ah  m  2  based  on  the  LiCoC^  electrode.  In  the 
base  case,  0.368  W  m  2-I<  is  used  for  h  to  represent  a  disadvanta¬ 
geous  heat-transfer  condition  [35],  Grid  resistance  ( RG )  [36]  is 
added  in  the  Dualfoil  program,  which  includes  the  total  resistance 
in  foils,  leads,  and  contacts  of  the  cell,  and  is  a  large  component  of 
internal  cell  resistance.  Typical  internal  resistance  is  in  the  order  of 
milliohms.  For  the  base-case  simulation,  RG  is  set  to  be 
0.8  x  1CT3  fi  m2,  and  RGext  is  set  to  be  0.2  x  1CT3  O  m2.  RGex t  is  the 
grid  resistance  outside  the  cell  sandwich  (but  not  Rioad)- 

Matrix  conductivity  is  assumed  to  be  high  value,  and  film 
resistance  reflects  the  resistance  of  SEI  layer  on  the  anode  [24],  The 
exchange  current  density  (io)  depends  on  the  local  instantaneous 
concentrations  as  computed  in  the  program. 

i0  =  Fkca“  (ct  -  cs)“cc““  [6] 


Table  4 

Activation  energy  values.3 


Activation  energy 

Value  ( Ea/R ,  units  of  K) 

Value  (J  mol  1) 

Solid-state  diffusion,  negative 

200 

1663 

Solid-state  diffusion,  positive 

900 

7483 

Electrolyte  conductivity  (k) 

1690 

13756 

Electrolyte  diffusion  (D) 

2000 

16628 

Negative  kinetics 

1800 

14965 

Positive  kinetics 

1800 

14965 

Film  resistance  of  anode 

-1800 

14965 

Film  resistance  of  cathode 

-1800 

14965 

a  Assumed  based  on  ref.  [48-52]. 
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In  some  simulation  literature,  only  initial  exchange  current 
density  values  are  given,  without  the  rate  constant  (fc).  Here  we  not 
only  show  k,  but  also  calculate  i0  value  for  the  readers'  convenience. 
The  values  of  reaction  rate  constants  are  set  to  be  relatively  high, 
and  the  two  insertion-electrode  reaction  kinetics  are  assumed  to  be 
facile,  and  thus  not  the  limiting  step.  The  exchange  current  varies 
with  concentration,  and  we  can  calculate  with  one  composition  as 
an  example.  When  c  =  1000  mol  m~3,  the  calculated  exchange 
current  density  (io)  at  298  I<  for  Lio.sCg  and  U0.6C0O2  by  Equation  (6) 
is  110  and  220  A  nrr2,  respectively.  Compared  with  values  used  in 
other  references,  these  values  are  relatively  high.  The  densities  in 
Table  2  are  crystalline  values  as  reported  in  the  literature  37], 

The  temperature  of  the  cell  can  be  calculated  by  an  energy 
balance  [19,20,38], 

-  J  E  'niUHjdv  -IV-  I2RGext  =  MCP§  +  h(T  -  Ta),  [7] 

V  l 


and  the  heat  produced  by  the  external  grid  resistance  (RGext)  has 
been  deducted.  L/h  is  the  enthalpy  potential  and  can  be  calculated 
by:  L/h  =  U-TdU/dT.  The  relations  of  U  (open-circuit  potential)  and 
dU/dT  with  lithium  concentration  for  UyCo02  and  MCMB  (LixCs) 
electrodes  [38]  are  displayed  in  the  Appendix. 

The  Dualfoil  program  has  been  developed  for  many  years  and 
various  electrolyte  systems  been  put  in  the  program  according  to 
our  research  [39,40].  For  the  state-of-art  lithium  ion  battery,  the 
most  common  electrolyte  system  is  alkyl  carbonate  based,  and  the 
Li-ion  battery  industry  does  not  have  as  many  choices  for  electro¬ 
lyte  as  there  are  for  positive  materials.  Table  3  lists  the  parameters 
for  the  electrolyte  used  in  our  simulation.  Here  the  commonly  used 
Bruggeman  expression  (e15)  is  used  to  correct  the  salt  diffusion 
coefficient  and  ionic  conductivity  for  the  base-case  simulation.  The 
lithium-ion  transference  number  and  thermodynamic  factor  are  set 
to  be  constant.  There  are  other  studies  on  electrolyte  [41—43], 

To  realize  fully  thermal-electrochemical  coupling,  the  depen¬ 
dence  of  physical  properties  on  temperature  must  be  considered. 
We  assume  an  Arrhenius-type  relationship 


>P  =  lL,oexP 


Ea(T  -  298)' 
298 RT 


[8] 


between  the  physical  properties  (Wo  represents  D,  k,  k,  Ds,  Rf  at 
298  K)  and  temperature.  Due  to  lack  of  abundance  of  direct 
experimental  data  [35,42,44—47],  these  activation-energy  values 
(Table  4)  are  assumed  based  on  classical  textbooks  [48—52], 


4.  Results  and  discussions 

4.2.  Superposition  integral  and  variable  Ds 

A  difficult  problem  in  modeling  research  existing  for  many  years 
is  how  to  simulate  a  variable  solid  diffusivity  (Ds)  not  only  based  on 
physics  and  experimental  property  measurements  but  also  in  a  short 
time.  Efforts  in  numerical  methods  have  been  summarized  by 
Ramadesigan  et  al.  [53]  and  Zeng  et  al.  [54].  Although  we  know  that 
Ds  varies  with  concentration  and  temperature  by  experiments,  an 
effective  mathematical  method  is  still  needed  to  treat  variable  Ds 
approximately  and  obtain  simulation  results  in  acceptable  accuracy. 

The  Dualfoil  model  has  two  modes  of  computation,  and  we  can 
set  mvdcl  or  mvdc3  as  1  in  the  input  file  to  run  the  program  with 
variable  Ds  (see  introduction  to  Dualfoil  on  our  website)  [14],  When 
Ds  does  not  vary  (over  time  and  position,  and  hence  over  temper¬ 
ature  and  composition),  it  is  appropriate  to  use  a  Duhamel  super¬ 
position  integral.  This  saves  greatly  in  computation  time.  If  Ds 


varies,  it  is  proper  to  use  a  finite-difference  treatment  in  the  solid 
particles.  This  is  much  slower,  and  it  can  run  into  convergence 
problems.  We  elected  to  use  the  faster  and  more  robust  super¬ 
position  integral  to  compute  variable  Ds,  and  in  practice  we  find  this 
method  gives  good  results.  This  approximate  approach  is  also 
adopted  by  other  researchers  to  get  rapid  simulation  [55], 

Fig.  1  gives  an  idea  of  the  error  incurred.  It  shows  current  and 
temperature  calculated  with  the  superposition  method  (curves  a 
and  b,  which  go  out  to  greater  capacity)  and  the  more  accurate 
finite-difference  method  (curve  c,  which  fails  earlier).  When  the 
activation  energies  for  Ds  are  zero,  there  is  less  actual  variation  of 
Ds,  and  the  superposition  method  should  become  comparable  in 
accuracy  to  the  finite-difference  method.  Fig.  1  compares  the  three 
numerical  modes.  The  original  Dualfoil4.2  model  using  the  Duha¬ 
mel  superposition  integral  is  strictly  valid  only  for  dealing  with 
constant  Ds  (case  a).  Dualfoil5.0  model  adds  a  pseudo-2D  procedure 
to  deal  with  variable  Ds  (case  c),  but  it  takes  longer  and  may  not 
converge.  Fig.  1  indicates  how  large  an  error  one  might  expect 
when  using  superposition  (case  b)  as  an  approximation  with  vari¬ 
able  Ds,  as  compared  with  case  c.  A  special  set  of  parameters  and 
constant  current  mode  are  adopted  to  get  as  much  data  for  case  c 
(long  simulation  time)  and  give  readers  visualized  comparisons 
between  the  three  numerical  modes.  To  get  faster  but  converged 
results,  in  the  following  simulations  we  use  the  superposition 
method  to  treat  with  variable  negative  and  positive  Ds  for  the 
modeling  purpose  of  very-high-current  discharge  mode  with  an 
appropriate  accuracy,  reducing  the  simulation  time  of  this  spherical 
diffusion  process  effectively  for  large-scale  simulations.  A  typical 
discharge  of  the  cell  under  constant  resistance  is  simulated  in  less 
than  20  min  using  this  approximation. 

4.2.  Base-case  simulation 

The  simulated  results  of  current  and  temperature-rise  behavior 
for  case  I  (base  case)  and  case  II  are  shown  in  Fig.  2(a)  and  (b).  The 
in-depth  analysis  of  the  two  cases  is  shown  in  Fig.  3(a)— (f).  Case  I  is 
the  base  case,  and  solution  diffusivity  (D)  is  dependent  on  tem¬ 
perature  with  an  activation  energy  ( Ea )  of  16.6  kj  mol-1  (as  shown 
in  Table  4).  For  comparison,  we  change  only  one  parameter  of  the 
base  case,  turning  off  the  thermal  activation  of  D,  which  is  case  II. 
Case  II  is  the  situation  that  solution  diffusivity  (D)  is  independent 
on  temperature,  with  zero  activation  energy,  maintaining  the  D 
value  at  298  K,  but  still  including  the  composition  dependence 
shown  in  Table  3. 


Capacity  (Ah/m2) 

Fig.  1.  Simulation  example  of  current  and  temperature-rise  behavior  versus  discharge 
capacity  for  comparison  of  the  three  numerical  modes:  a.  Superposition  method  for 
constant  Ds;  Ea/R  for  negative  and  positive  Ds  is  0.  b.  Superposition  method  for  variable 
negative  and  positive  Ds;  Ea/R  for  negative  and  positive  Ds  is  200  and  900  K,  respec¬ 
tively.  c.  Finite-difference  method  for  variable  negative  and  positive  Ds;  Ea/R  for 
negative  and  positive  Ds  is  200  and  900  K,  respectively. 
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Utilization  of  LiyCoOz 

Fig.  2.  (a)  Comparison  of  current  and  temperature  simulation  results  versus  time  for 
case  I  and  case  II.  (b)  Comparison  of  cell  potential  and  temperature  simulation  results 
versus  utilization  of  positive  material  [calculated  by  (y-0.6)/0.4].  Case  I  is  the  base  case, 
and  solution  diffusivity  (D)  is  dependent  on  temperature  with  an  activation  energy  ( Ea ) 
of  16.6  kj  mol-1.  Case  II  is  the  situation  that  D  is  independent  of  temperature,  keeping 
the  value  at  298  K.  Similar  information  is  plotted  in  Figure  a  and  b  to  show  how  it 
appears  with  different  abscissas. 


The  simulated  results  in  Fig.  2(a)  and  (b)  for  the  two  cases  are 
very  different.  This  is  due  to  the  difference  of  thermal- 
electrochemical  coupling  effect  and  the  difference  of  solution 
diffusion  coefficient  (D).  For  both  cases,  the  zero-time  peak  cur¬ 
rent  is  controlled  by  the  sum  of  the  resistance  of  the  external  load 
and  the  effective  resistance  of  the  electrochemical  cell.  The  cur¬ 
rent  distribution  within  the  porous  electrodes  is  controlled  by  the 
electrode  kinetics  and  the  conductivities  of  the  solution  and  solid 
phases  [56,57],  For  times  larger  than  zero,  the  current  initially 
declines  rapidly  due  to  mass-transport  limitations  of  lithium  ions 
in  solution  and  the  rise  in  surface  concentration  of  lithium  on  the 
positive  electrode  surface  (see  Fig.  3).  The  negative  electrode  is 
not  showing  significant  limiting  behavior  characteristics  since 
the  lithium-ion  concentration  is  rising  within  the  porous 
electrode. 

For  the  base  case,  as  temperature  rises,  the  solution  and  solid- 
state  diffusion  coefficients  begin  to  rise,  resulting  in  improved 
mass  transport  of  lithium  in  the  solution  and  solid  phases  leading  to 
a  reversal  of  the  decline  in  current.  The  current  continues  to  in¬ 
crease  as  mass  transport  increases.  During  this  period,  the  lithium 
ion  concentration  in  the  solution  increases  while  the  surface  con¬ 
centration  of  lithium  begins  to  approach  unity  along  the  length  of 
the  pore.  When  the  peak  in  the  current  is  reached,  mass  transport 
of  lithium  into  the  solid  particle  becomes  limiting  (the  surface 
concentration  is  beginning  to  approach  the  maximum  concentra¬ 
tion  of  particle),  and  the  current  begins  a  continuous  decline  over  a 
period  of  minutes.  The  current  plot  shows  a  sharp  peak  after  the 
quick  drop,  and  the  cutoff  potential  (0.0001  V)  for  simulation  is 


reached  rapidly.  The  current  has  its  maximum  value  of  2514  A  nT2 
at  the  beginning  of  discharge  (point  b,  imax).  drops  sharply  to  the 
valley  value  of  823.3  A  itT2  (point  c,  iVaiiey).  climbs  to  the  sharp  peak 
value  of  1932  AnT2  (point  e,  ipeak).  and  then  decreases  gradually  to 
near  zero.  However,  for  case  II,  due  to  no  thermal  activation  of  the 
solution  diffusion  coefficient,  the  solution  limitation  has  no  chance 
to  decrease,  and  the  second  sharp  peak  disappears;  the  current 
shows  a  gentle  slope  after  the  rapid  drop  to  300  A  nT2,  and  the 
cutoff  potential  is  reached  in  a  longer  time.  The  cell  temperature 
rises  more  rapidly  for  the  base  case,  reaching  120  °C  in  12  s,  while 
the  cell  temperature  reaches  120  °C  in  27  s  for  case  II. 

As  shown  in  Fig.  2(b),  for  the  same  utilization,  the  base  case 
shows  a  lower  cell  temperature  rise  than  case  II  due  to  its  lower 
solution  mass-transport  limitation;  for  the  same  discharge  time, 
the  base  case  shows  a  higher  utilization  (a  faster  discharge).  For  the 
base  case,  most  of  the  capacity  is  given  out  at  the  sharp  peak;  the 
average  utilization  of  positive  material  reaches  0.6  at  30.5  s.  For 
case  II,  most  capacity  is  given  out  at  the  gentle  slope;  the  average 
utilization  of  positive  material  reaches  only  0.23  at  30.5  s.  The 
temperature-rise  rate  is  approximately  proportional  to  the 
discharge  capacity  (or  utilization  of  the  positive  material)  for  both 
cases.  The  reason  for  higher  maximum  temperature  for  case  II  is 
that  the  utilization  is  higher. 

Heat  can  be  produced  by  electrochemical  reaction,  chemical  side 
reactions,  and  Joule  heating  inside  the  battery.  In  these  simulations, 
chemical  side  reactions  are  not  considered,  and  only  the  normal 
electrochemical  reaction  is  included.  At  the  starting  moment  of 
discharge,  the  internal  cell  resistance  is  composed  of  RG 
(0.8  x  10~3  Q  m2),  other  solid  and  solution  ohmic  resistance,  and 
electrochemical-kinetic  resistance;  the  load  resistance  is 
0.2  x  10~3  Q  m2;  Uocp  is  about  3.8  V;  and  thus  the  current  density  is 
very  high  (more  than  2500  A  nT2),  and  Joule  heat  (Q  =  i2Rt)  is 
dominant.  The  cell  temperature  rises  quickly  to  more  than  50  °C  at 
4.2  s  for  both  cases.  The  quick  drop  of  current  for  both  cases  results 
from  the  solution  mass-transport  limitation;  the  high  current 
certainly  depletes  the  solution  in  the  back  of  the  positive,  and  the 
surface  of  the  active  material  in  the  front  of  the  positive  becomes 
nearly  saturated.  The  solution  mass-transport  resistance  is 
dominant. 

To  understand  better  the  difference  between  base  case  and  case 
II,  the  solution  concentration  profiles  are  simulated  and  shown  in 
Fig.  3(a)  and  (b).  For  the  base  case,  thanks  to  the  thermal  activation 
of  D,  the  solution  gradients  decrease  gradually,  the  solution  then 
becomes  more  and  more  uniform,  and  the  solution  limitation  dis¬ 
appears  finally  at  point  f.  That  is  the  reason  for  the  sharp  current 
peak  appearing.  But  for  case  II,  without  thermal  activation  of  D,  the 
solution  mass  transfer  has  no  chance  to  increase.  The  solution 
gradients  stay  large,  and  the  solution  limitation  exists  all  the  time. 
That's  why  the  current  declines  at  a  gentle  slope.  When  solution 
limitations  appear,  the  solution  mass-transport  resistance  will  be 
large,  which  can  account  for  the  quick  drop  at  4.20  s  for  both  cases. 

As  shown  in  Fig.  3(c)  and  (d),  for  both  cases  at  4.2  s,  the  solid 
surface  concentration  in  the  front  of  the  positive  electrode  reaches 
near  saturation  (y  approaches  1)  due  to  high  electrolyte  concen¬ 
tration  and  favorable  potential  near  the  separator;  the  solid  surface 
concentration  in  the  back  of  positive  electrode  is  much  smaller  due 
to  the  low  electrolyte  concentration  and  unfavorable  potential.  As 
the  solid  surface  concentration  nears  saturation,  the  solid  diffu¬ 
sivity  is  also  a  limitation  for  both  cases.  For  example,  as  for  the  base 
case  at  30.5  s,  the  solution  concentration  is  uniform  (solution-mass 
transport  limitation  disappears)  due  to  the  thermal- 
electrochemical  coupling  effect;  the  y  value  in  the  LiyCo02  surface 
is  more  than  0.98  (nearly  saturated),  but  the  average  y  value 
through  the  whole  positive  electrode  is  only  0.8362,  which  means 
that  the  core  of  the  solid  particle  still  has  a  lot  of  vacancies  for 
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Fig.  3.  (a)  and  (b)  Simulated  solution  concentration  profiles  for  case  I  and  case  II  in  different  depths  of  the  cell,  (c)  and  (d)  Simulated  solid  surface-concentration  profiles  for  case  I 
and  case  II  in  different  depths  of  the  LiyCo02  positive  electrode,  (e)  and  (f)  Simulated  solution-solid  interface  transfer  current  profiles  for  case  1  and  case  II  in  different  depths  of  the 
LiyCo02  positive  electrode. 


lithium  ions,  and  the  discharge  is  controlled  by  lithium  ions 
diffusing  from  the  surface  into  the  core. 

According  to  the  Butler— Volmer  equation,  transfer  current  j  is 
determined  by  solution  concentration,  solid  surface  concentration, 
and  overpotential.  Fig.  3(e)  and  (f)  shows  the  distribution  of 
transfer  current  across  the  positive  electrode.  For  both  cases,  at  the 
front  of  the  positive  electrode,  a  solution  limitation  does  not  exist, 
the  solid  surface  approaches  saturation,  and  transfer  current  shows 
a  constant  value,  indicating  a  steady  discharge  state  controlled  by 
solid  diffusion.  For  the  base  case,  due  to  the  temperature  rise  and 
increase  of  D  and  Ds,  the  electrolyte  diffuses  to  the  back  of  the 
positive  electrode  gradually,  and  the  transfer  current  increases.  At 
30.5  s,  the  transfer  current  is  uniform  and  almost  constant.  For  case 
II,  at  30.5  s,  the  transfer  current  at  the  back  of  the  positive  electrode 
is  still  very  small.  This  is  because  the  solution  concentration  is  still 
near  zero,  lacking  the  thermal  activation  of  D.  But  for  case  II,  from 
point  c  to  f,  the  transfer  current  shows  the  trend  of  becoming  more 


uniform;  this  must  be  due  to  the  thermal  activation  of  solid  diffu¬ 
sion,  increasing  Ds. 

From  the  above  discussion,  it  is  found  that  the  temperature 
dependence  of  D  has  a  large  impact  on  the  discharge  behavior. 
Therefore,  we  calculated  the  diffusion  coefficient  by  Equation  (8). 
For  the  base  case,  from  point  c  to  e,  the  temperature  varies  from 
65.5  to  185.38  °C;  Ds,p  is  doubled,  and  D  is  increased  to  4.5  times. 

4.3.  Simulation  of  sensitivity  tests 

In  the  following,  a  series  of  sensitivity  tests  is  conducted  to 
investigate  the  effect  of  each  parameter  on  the  cell  temperature. 
Without  doing  expensive  and  dangerous  battery  safety  tests  [58], 
with  the  help  of  numerical  simulations,  we  can  explore  the  cell 
performance  in  a  wide  range,  obtain  deep  understanding  of  over¬ 
heating  inside  the  cell,  and  interpret  many  experimental 
phenomena. 
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Mass-transport  resistance  is  a  large  source  of  internal  resistance 
determining  the  heat  production.  Firstly,  we  change  parameters 
related  to  solution  mass  transport.  Secondly,  we  change  parameters 
related  to  solid  mass  transport.  Then,  we  change  parameters  related 
to  heat  transfer  to  the  environment. 

The  main  feature  of  the  lithium  ion  battery  is  its  porous  struc¬ 
ture.  Electrolyte  diffuses  inside  the  pores  (a  limited  space)  formed 
by  solid  material.  The  effective  values  of  transport  properties  (Deff 
and  Keff)  in  porous  electrodes  are  smaller  than  the  intrinsic  value 
which  is  measured  in  a  free  solution.  The  Bruggeman  relation 
(Equation  (9))  is  often  used  to  calculate  the  effective  value. 

Deff  =  De“  [9] 

The  porous  networks  of  negative  electrode,  positive  electrode, 
and  separator  are  different.  In  the  Dualfoil  model,  the  Bruggeman 
exponent  (a)  is  now  set  separately  for  the  negative  electrode 
(Exbrugl),  the  separator  (Exbrug2),  and  the  positive  electrode 
(Exbrug3).  The  Bruggeman  exponent  has  been  shown  to  range  from 
1.5  to  3  from  the  information  in  the  literature  and  from  other 
battery  systems  [59],  In  the  base  case, 
Exbrugl  =  Exbrug2  =  Exbrug3  =  1.5.  (In  reference  12,  a  factor  of  e  is 
taken  out  already  for  D,  and  the  exponent  would  then  be  0.5  instead 
of  1.5.  In  Dualfoil,  a  definition  like  Equation  (9)  is  used  for  Deff  as 
well  as  Keff.) 

Fig.  4(a)  indicates  that  decreasing  the  positive  electrode 
porosity  increases  the  valley  current,  and  the  sharp  peak  ap¬ 
pears  later.  Fig.  4(b)  indicates  that  increasing  Exbrug3  decreases 
the  valley  current,  and  the  sharp  peak  appears  later  with  a 
smaller  peak  current.  For  a  real  battery,  the  porosity  can  vary 
over  only  a  narrow  range,  but  the  Bruggeman  exponent  can  be 


Time  (s) 

Fig.  4.  (  a)  Simulated  current  and  temperature  results  for  different  porosities  of  positive 
electrode,  (b)  Simulated  current  and  temperature  results  for  different  Bruggeman 
exponents  of  positive  electrode  (Exbrug3). 


varied  over  a  wide  range  due  to  material  morphology  [59], 
Compared  with  the  porosity,  the  Bruggeman  exponent  has  a 
larger  effect  on  a  and  D. 

The  solid-diffusion  process  inside  the  particles  is  determined 
by  Ds  p  and  rp  Many  parameters  such  as  doping  element,  lattice 
defects,  lithium  concentration,  particle  morphology,  and  elec¬ 
trode  microstructure  can  affect  the  effective  solid-diffusion  co¬ 
efficient  of  LiyCo02  particles  in  the  porous  electrode.  For 
simplicity,  in  this  simulation,  Ds,p  at  25  °C  is  assumed  to  be  the 
effective  value  and  is  adjustable.  The  diffusion  distance  for 
lithium  ions  inside  the  positive  spherical  particle  is  the  radius 
(rp),  which  can  be  adjustable  in  a  wide  range.  Fig.  5(a)  shows  that 
decreasing  Ds  p  low  enough  causes  the  current  peak  to  disappear; 
the  temperature  rising  rate  and  the  maximum  temperature  is 
decreased  substantially.  Fig.  5(b)  shows  that  increasing  the  radius 
of  LiyCo02  particles  causes  the  current  peak  and  the  cell  tem¬ 
perature  to  decrease  gradually.  If  the  LiyCo02  particle  is  big 
enough,  the  current  peak  will  disappear  too,  and  the  cell  tem¬ 
perature  will  be  below  120  °C.  Fig.  5(c)  shows  the  current  peak  of 
a  series  of  sensitivity  tests  on  combinations  of  Ds  p  and  rp.  (The 
temperature  rise  behavior  is  not  shown  here.)  The  red  line  is  the 
boundary  for  peak  and  no  peak.  For  very  large  particles  and  very 
low  Ds,p,  the  LiyCo02  surface  should  be  saturated  rapidly  during 
the  high-current  discharge,  and  then  the  discharge  is  controlled 
by  solid  diffusion  with  a  small  transfer  current.  There  should  be 
no  current  peak,  and  the  cell  temperature  will  not  rise  to  more 
than  120  °C,  and  no  thermal  runaway  should  happen.  The 
behavior  is  similar  to  the  effect  of  D,  where  shutting  off  the 
thermal  activation  of  solution  diffusivity  causes  the  peak  to 
disappear  (see  Fig.  2).  Through  extensive  simulations,  we  can 
understand  many  experimental  results  for  short-circuit  tests  and 
that  some  high-power  type  Li-ion  batteries,  which  commonly 
have  smaller  particles,  show  a  higher  possibility  of  overheating 
during  nail  penetration  tests  [5—7], 

Under  the  same  heat-transfer  conditions,  the  temperature-rise 
rate  for  current-peak  appearance  is  higher  than  that  for  no 
current-peak  appearance.  By  choosing  the  proper  electrode  mate¬ 
rial  and  battery  design,  the  cell  temperature-rise  rate  under  very 
high  current  discharge  can  be  decreased  a  little. 

A  cell  consisting  of  only  one  positive,  negative,  and  separator 
and  with  external  turbulent  air  flow  might  correspond  to  a  heat- 
transfer  coefficient  of  36.8  W  nT2  K.  A  cell  with  a  few  more  elec¬ 
trodes  and  free  convection  on  the  outside  might  correspond  to 
3.68  W  nrT2  K.  For  a  larger,  perhaps  spirally  wound  cell,  a  value  of 
0.368  W  uT2  K  might  give  a  better  representation  of  the  average 
temperature,  which  can  be  substantially  higher  toward  the  core  of 
the  cell.  With  higher  h,  heat  can  go  out  of  the  cell  faster;  then  the 
temperature-rise  rate  decreases,  the  maximum  cell  temperature 
decreases,  and  the  peak  current  recovered  by  thermal  activation 
will  decrease,  and  even  disappear. 

Note  that  the  maximum  cell  temperature  is  only  94.6  °C  for  h  of 
36.8  W  nrT2  I<  ',  a  very  small  value  compared  with  that  for  small  h. 
Because  the  solution  and  solid-diffusion  limitations  always  exist 
under  very-high-current  discharge,  the  heat  produced  by  internal 
resistance  (mainly  mass-transport  resistance)  is  large.  If  the  heat 
cannot  be  dissipated  rapidly  (small  h),  the  thermal-electrochemical 
coupling  effect  will  speed  up  the  discharge  process;  the  full  ca¬ 
pacity  will  be  released  rapidly  and  cause  the  cell  temperature  to 
reach  120  °C  more  rapidly  (See  Fig.  6). 

The  simulation  results  for  different  h  values  tell  us  that  it  is 
effective  to  decrease  the  cell  temperature  by  increasing  the  heat- 
transfer  coefficient.  For  a  large  and  thick  lithium-ion  battery,  the 
outer  side  of  the  cell  has  a  big  h,  but  the  center  of  the  cell  has  a  small 
h.  Heat  is  hard  to  dissipate  from  the  center  of  the  cell.  For  a  lithium- 
ion  battery  used  in  electric  vehicles,  forced-convection  cooling 
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Fig.  5.  (a)  Simulated  current  and  temperature  results  for  different  Ds>p  (b)  Simulated  current  and  temperature  results  for  different  radius  of  LiyCo02  particles,  (c)  Simulated  relative 
current-peak  areas  of  different  combinations  of  Ds  p  and  rp.  (Solid  circles'  area  shows  relative  size  of  corresponding  peak  area  qualitatively,  and  open  circles  mean  no  peak  for 
current.  The  dashed  line  is  the  boundary  of  peak  and  no  peak.) 


Utilization  of  Li  CoO„ 

y  2 

Fig.  6.  Simulated  current  and  temperature-rise  behavior  for  different  heat-transfer 
conditions. 


conditions  should  be  very  helpful  to  decrease  the  surface  temper¬ 
ature  of  the  cell.  But  we  still  need  proper  battery  design  to  dissipate 
heat  from  the  center  of  the  cell  effectively. 

4.4.  Shut-down  separator  and  side  reactions 

To  prevent  the  cell  temperature  from  reaching  elevated  values 
observed  in  commercial  cells  and  shown  in  Fig.  7,  a  “shut-down” 
separator  is  used  in  which  the  pores  begin  to  close  when  the 
temperature  approaches  approximately  127  °C.  The  commonly 
used  separator  is  composed  of  a  porous  poly(ethylene)  (PE)- 
polypropylene  (PP)  bilayer  or  a  PP-PE-PP  trilayer  structure,  which 
can  shut  itself  down  by  closing  the  pores  when  the  cell  temper¬ 
ature  is  higher  than  the  melting  point  of  PE  (about  127  °C)  [60], 
After  the  shut-down,  lithium-ions  cannot  pass  through  the  pores 
to  reach  the  cathode,  and  the  normal  electrochemical  discharge 
process  stops.  However,  should  the  heating  rate  become  large 
[60  ,  closure  of  the  pores  may  not  be  fast  enough  and  the  tem¬ 
perature  may  increase  beyond  the  shut-down  rated  value 
(measured  at  1  °C  s_1)  before  the  pores  are  sufficiently  closed  to 
terminate  the  flow  of  current.  Fig.  7  shows  that  the  simulated 
behavior  of  the  cell  has  an  associated  heating  rate  >60  °C  s.  This 
condition  suggests  that  the  temperature  can  easily  reach  tem¬ 
peratures  exceeding  127  °C.  This  condition  will  result  in  acceler¬ 
ated  heat  generation  caused  by  elevated  rates  of  side  reactions 
[61],  Under  this  condition  the  possibility  of  a  thermal  runaway 
can  occur. 

Side  reactions  are  increased  by  rapidly  rising  temperatures.  If 
heat  generation  exceeds  that  can  be  removed  from  the  cell,  the 
temperature  will  continue  to  rise.  The  lithium  ion  cell  can 
perform  well  under  moderate  temperature  and  current  density 
only  because  the  kinetics  of  any  side  reaction  is  slow  due  to  the 
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Time  (s) 

Fig.  7.  Simulated  current  and  temperature  curves  versus  time  for  R|0,l(l  =  0.0002  Q  m2 
with  and  without  shut-down  separator,  h  =  0.368  W  m  2  K  '. 


formation  of  a  solid-electrolyte  interphase  (SEI)  [62—67].  The  side 
reactions  of  Li-ion  cells  reported  in  the  literature  include  the  SEI 
film  growing  and  decomposing  on  carbon  negative  material, 
solvent  reduction,  and  decomposition  of  positive  material,  salt 
decomposition,  and  fracture  of  positive  and  negative  solid  ma¬ 
terial;  continued  rise  in  temperature  results  in  a  cell-pressure 
increase  from  a  rise  in  vapor  pressure  and  gas  released  by  side 
reactions  [68—70],  If  the  pressure  becomes  too  high,  a  pressure- 
relief  valve  will  be  activated  [71],  Even  if  the  separator  is  shut 
down,  and  the  normal  discharge  is  stopped,  some  side  reactions 
still  can  happen  in  the  negative  or  positive  electrode  domains 
[72,73],  More  novel  ideas  such  as  thermo-responsive  polymer 
added  in  the  anode  to  aid  shutdown  are  needed  to  improve 
battery  safety  74],  The  simulation  of  combining  the  shutdown 
separator  and  the  side  reactions  for  lithium-ion  cell  by  the 
Dualfoil  model  is  being  studied  now  and  will  be  reported  in  the 
future. 

5.  Conclusions 

With  a  radial  mass-transport  coefficient  to  modify  the  thermal- 
electrochemical-  coupled  Dualfoil  model,  simulation  of  very-high- 
current  discharge  of  a  lithium-ion  cell  is  achieved.  This  modifica¬ 
tion  not  only  allows  the  program  to  run  at  virtually  all  currents, 
but  it  also  dramatically  reduces  the  running  time,  allowing  for 
timely  exploration  of  the  impact  of  changes  in  the  various  cell 
parameters.  For  small-resistance-load  discharge,  the  current  is 
high,  and  the  cell  temperature  rise  is  large.  The  solution  concen¬ 
tration  in  the  back  of  the  positive  electrode  becomes  quite  small, 
and  the  front  of  the  positive  becomes  saturated  with  lithium  at  the 
surface  of  the  particles.  The  current  drops  because  it  is  now  limited 
by  diffusion  within  the  particles  at  the  front  of  the  electrode. 
Under  certain  circumstances,  this  drop  in  cell  potential  and  current 
occurs  even  though  there  is  a  lot  of  capacity  in  the  back  of  the 
positive  and  some  remaining  at  the  center  of  the  particles  near  the 
separator.  Meanwhile,  the  large  rise  in  cell  temperature  leads  to  an 
increase  in  the  diffusion  coefficients  both  in  the  solution  and  in  the 
particles,  and  there  is  a  dramatic  increase  in  cell  potential  and 
current,  permitting  the  rest  of  the  cell  capacity  to  be  discharged 
rapidly.  Sensitivity  tests  indicate  that  Li-ion  cells  composed  of 
small  particles  with  large  effective  solid  diffusivity  show  more 
rapid  temperature  rise  under  high-current  discharge.  Special  bat¬ 
tery  design  is  needed  to  achieve  effective  heat  transfer  outside  the 
cell  center  and  control  the  cell  temperature  in  a  proper  range.  An 
effective  mechanism  needs  to  be  developed  to  shut  down  the 
exothermic  side  reaction  inside  the  cell  completely  when  the  cell 
is  overheating. 
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Appendix 


The  equations  used  in  the  Dualfoil  model  are  listed  in  Table  5. 
The  open-circuit  potential  dU/dT  profiles  used  in  the  simulations 
are  based  on  the  results  from  Refs.  [38,75,76]  modified  slightly,  and 
displayed  in  Fig.  8.  The  fitting  expressions  are  presented  in  Equation 

(10)— (13). 

For  the  positive  LiyCo02  electrode. 


Uy  = 


2.16216  +  0.07645  tanh(30.834  -  54.4806y) 
+  2.1581  tanh(52.294  -  50.294y) 

-  0.14169  tanh(l  1.0923  -  19.8543y) 

+  0.2051  tanh(1.4684  -  5.4888y) 


+  0.2531  tanh 
-  0.02167  tanh 


-y  +  0.56478\ 
0.1316  ) 

y-0.525\ 
0.006  ) 


[10] 


d Uy  -0.19952  +  0.92837y  -  1.36455y2+  0.61 154y3 

~dT~l  -  5.66148y  +  11.47636y2-  9.82431y3  +  3.04876y4 

[11] 

For  the  negative  LixC6  (MCMB)  electrode. 

Ux  =  0.194  +  i.5e(-120'0x>  +  0.0351  tanh* “  0  286 


0.083 
x-  0.9233 


v  _  n  849 

~  0-0045  tanh^^- 0.035  tanh  Q  05 
-  0.0147  tanh*  O'.5  -  0.102  tanh* -  — 94 


0.034 


0.142 
x-  0.124 


-  0  022  tanh*^- O.OH  tanh 


+  0.0155  tanh 


x-  0.105 
0.029 


dUx  _  A 
d T~B 


[12] 

[13] 


Table  5 

Equations  used  in  the  Dualfoil  model. 


Equation  description  Equation 


Electrolyte  material 
balance 


«2  f  =V-(e2DVc) 


fa(1  ij) 


Intercalant  material 
balance 
Liquid-phase 
Ohm's  law 
Solid-phase 
Ohm's  law 
Butler— Volmer 
kinetics 

Exchange  current 
density 

Charge  conservation 


acs  _  i  e  n  r2  dcs  \ 

W~Wd r  [Us'  W J 

i2  =  -/cV<*> 2+^(1  -  t$)Vln(f±c) 

l’l  =  —  crVdh 

1  =jT  =  ‘o  [exp (jffvs)  -exp^-^7,s)j 
i0  =  Fkc“°(ct-cs)“‘cZ‘ 


V  i2  =  aFj 


j  denominator 
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Fig.  8.  The  open-circuit  potential  and  dil/dT  profiles  for  MCMB  (LixC6)  and  LiyCo02  electrode.[38,75,76]. 


A  =  0.00527  +  3.29927x  -  91.79326x2  +  1004.91101x3  -  5812.27813x4  +  19329.75490x5  -  37147.89470x6 
+  38379. 18127x7  -  16515.05308x8 

B  =  1  -  48.09287X  +  1017.23480x2  -  1 0481. 8041 9x3  +  59431. 30001x4  -  195881.64880x5  +  374577.31520x6 
-  385821. 16070x7  +  165705.85970x8 


List  of  symbols 

him 

limiting  transfer  current,  A  m~2 

in, I 

transfer  current  at  the  node  (n,  /),  A  itT2 

a 

active  interfacial  area  per  unit  electrode  volume,  nr1 

j 

transfer  flux,  mol  m_2.s 

Cp 

heat  capacity,  J  kg-1  1< 

k 

insertion  reaction  constant,  m2  5  mol-0'5  s_1 

c 

electrolyte  concentration,  mol  m  3 

M 

mass  density  of  cell,  kg  itT2 

co 

electrolyte  concentration  at  the  electrode  surface, 

R 

gas  constant,  8.314  J  mol-1  I< 

mol  m~3 

r 

radius  of  the  pore,  m 

cs 

lithium  ion  concentration  in  solid,  mol  nrT3 

rP 

radius  of  the  positive  electrode,  m 

ct 

maximum  lithium  ion  concentration  in  solid,  mol  nT3 

Rf 

film  resistance  fi  m2 

D 

bulk  electrolyte  diffusion  coefficient,  m2  s_1 

T 

cell  temperature 

Deff 

effective  solution  diffusion  coefficient,  m2  s_1. 

ra 

ambient  air  temperature,  1< 

Ds 

solid  diffusion  coefficient,  m2  s_1 

u 

open-circuit  potential,  V 

Ds,p 

solid  diffusion  coefficient  of  positive  electrode,  m2  s_1 

UH 

enthalpy  potential,  V 

Ea 

activation  energy,  J  mol-1 

V 

volume  of  electrode 

F 

Faraday's  constant,  97485  C  mol-1 

Initial  stoichiometric  value  for  Li*C6 

f± 

salt  activity  coefficient 

y0 

Initial  stoichiometric  value  for  LiyCoC>2 

h 

heat  transfer  coefficient,  W  m  2  1< 

a 

Bruggeman  exponent 

1 

cell-sandwich  current  density,  A  uT2 

aa 

anodic  transfer  coefficient 

i 

transfer  current,  A  ur2 

ac 

cathodic  transfer  coefficient 

h 

solid-phase  current  density,  A  m~2 

8 

electrode  porosity  or  volume  fraction  of  solution 

h 

solution-phase  current  density,  A  m~2 

€f 

volume  fraction  of  inert  filler 

b 

exchange  current  density,  A  m  2 

K 

ionic  conductivity,  S/m 
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ris  surface  over-potential,  V 

a  electronic  conductivity,  S/m 

W  temperature-dependent  physical  properties 

W o  temperature-dependent  physical  properties  at  298  K 
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